# Replication File for "Combining List Experiment and Direct Question Estimates of Sensitive Behavior Prevalence"
# Peter M. Aronow, Alexander Coppock, Forrest W. Crawford, Donald P. Green
# Placebo Test I Power Calculations

rm(list=ls())
# Set your working directory to ACCG_list_replication_archive
setwd()

### Bring in functions
source("programs/ACCG_list_source.R")
library(lattice)
library(parallel)

### Prepare for Parallel
cl <- makeCluster(getOption("cl.cores", 4))
clusterEvalQ(cl, source("programs/ACCG_list_source.R")) 

#parameters for all figures
sims <- 1000
col.l <- colorRampPalette(c('white', 'black'))(100)
col.l[c(75:85)] <- c("red4","red3","red3","red2", "red2", "red","red2","red2","red3", "red3","red4")
depth.breaks <- do.breaks(c(0,1), 100)

# Figure A1: Power of Placebo Test I to Detect Violations of Assumptions 1 - 3

# Panel A1-1: Violated assumption: no false confessions
admitter_Ns <- seq(100,1000,by=50)
falseconfessor_props <- seq(0,1,by=.01)

powerz <- NULL
admitter_Nz <- NULL
falseconfessor_propz <- NULL

for(i in 1:length(admitter_Ns)){
  for(j in 1:length(falseconfessor_props)){
    N.trueadmitter <- floor(admitter_Ns[i]*(1-falseconfessor_props[j]))
    N.falseconfessor <- admitter_Ns[i] - N.trueadmitter
    N.withholder <- 200       #Engage = Yes, Say So = No
    N.innocent <- 200         #Engage = No, Say So = No
    clusterExport(cl, c("N.trueadmitter","N.falseconfessor", "N.withholder","N.innocent" ))
    power <- mean(parSapply(1:sims, cl=cl, FUN=function(i,...) { x <-  listsimulator(N.trueadmitter=N.trueadmitter,N.falseconfessor=N.falseconfessor,N.withholder=N.withholder,N.innocent=N.innocent,baseline_binomial_prob=.4)$reject_null}))
    admitter_Nz <- c(admitter_Nz, admitter_Ns[i])
    falseconfessor_propz<- c(falseconfessor_propz, falseconfessor_props[j])
    powerz <- c(powerz, power)
  }
}

graph.frame.1 <- data.frame(x=falseconfessor_propz, y=admitter_Nz, z = powerz)

powerfig1 <- levelplot(z~x*y, graph.frame.1, cuts=20, col.regions=col.l,
                   colorkey=FALSE,
                   at=depth.breaks,
                   ylab = "Total Number of Admitters",
                   xlab = "Proportion of false confessors",
                   main = "Violated assumption: no false confessions",
                   legend = 
                     list(right = 
                            list(fun = draw.colorkey,
                                 args = list(key = list(col = col.l, at = depth.breaks),
                                             draw = FALSE)))
)

pdf("placebopowerfig1.pdf")
powerfig1
dev.off()

# Panel A1-2: Violated assumption: no liars
admitter_Ns <- seq(100,1000,by=50)
listliar_props <- seq(0,1,by=.01)

powerz <- NULL
admitter_Nz <- NULL
listliar_propz <- NULL

for(i in 1:length(admitter_Ns)){
  for(j in 1:length(listliar_props)){
    N.trueadmitter <- admitter_Ns[i]
    N.falseconfessor <- 0
    N.withholder <- 200       #Engage = Yes, Say So = No
    N.innocent <- 200         #Engage = No, Say So = No
    listliar_prop <- listliar_props[j]
    clusterExport(cl, c("N.trueadmitter","N.falseconfessor", "N.withholder","N.innocent","listliar_prop" ))
    power <- mean(parSapply(1:sims, cl=cl, FUN=function(i,...) { x <-  listsimulator(N.trueadmitter=N.trueadmitter,N.falseconfessor=N.falseconfessor,N.withholder=N.withholder,N.innocent=N.innocent,baseline_binomial_prob=.4,pr.listliar=listliar_prop)$reject_null}))
    admitter_Nz <- c(admitter_Nz, admitter_Ns[i])
    listliar_propz<- c(listliar_propz, listliar_props[j])
    powerz <- c(powerz, power)
  }
}

graph.frame.2 <- data.frame(x=listliar_propz, y=admitter_Nz, z = powerz)
powerfig2 <- levelplot(z~x*y, graph.frame.2, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Total Number of Admitters",
                       xlab = "Proportion of liars on the treatment list",
                       main = "Violated assumption: no liars",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

pdf("placebopowerfig2.pdf")
powerfig2
dev.off()

# Panel A1-3: Violated assumption: no design effects
admitter_Ns <- seq(100,1000,by=50)
design_affected_props <- seq(0,1,by=.01)

powerz <- NULL
admitter_Nz <- NULL
design_affected_propz <- NULL

for(i in 1:length(admitter_Ns)){
  for(j in 1:length(design_affected_props)){
    N.trueadmitter <- admitter_Ns[i]
    N.falseconfessor <- 0
    N.withholder <- 200       #Engage = Yes, Say So = No
    N.innocent <- 200         #Engage = No, Say So = No
    design_affected_prop<- design_affected_props[j]
    clusterExport(cl, c("N.trueadmitter","N.falseconfessor", "N.withholder","N.innocent","design_affected_prop" ))
    power <- mean(parSapply(1:sims, cl=cl, FUN=function(i,...) { x <-  listsimulator(N.trueadmitter=N.trueadmitter,N.falseconfessor=N.falseconfessor,N.withholder=N.withholder,N.innocent=N.innocent,baseline_binomial_prob=.4,pr.designaffected=design_affected_prop)$reject_null}))
    admitter_Nz <- c(admitter_Nz, admitter_Ns[i])
    design_affected_propz<- c(design_affected_propz, design_affected_props[j])
    powerz <- c(powerz, power)
  }
}

graph.frame.3 <- data.frame(x=design_affected_propz, y=admitter_Nz, z = powerz)

powerfig3 <- levelplot(z~x*y, graph.frame.3, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Total Number of Admitters",
                       xlab = 'Proportion of "Design Affected"',
                       main = "Violated assumption: no design effects",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

pdf("placebopowerfig3.pdf")
powerfig3
dev.off()

# Panel A1-4: Variability of responses to non-sensitive items
admitter_Ns <- seq(100,1000,by=50)
baseline_probs <- seq(0,1,by=.01)

powerz <- NULL
admitter_Nz <- NULL
baseline_probz <- NULL

for(i in 1:length(admitter_Ns)){
  for(j in 1:length(baseline_probs)){
    N.trueadmitter <- floor(admitter_Ns[i]*(1-.20))
    N.falseconfessor <- admitter_Ns[i] - N.trueadmitter
    N.withholder <- 200       #Engage = Yes, Say So = No
    N.innocent <- 200         #Engage = No, Say So = No
    baseline_prob <- baseline_probs[j]
    clusterExport(cl, c("N.trueadmitter","N.falseconfessor", "N.withholder","N.innocent","baseline_prob" ))
    power <- mean(parSapply(1:sims, cl=cl, FUN=function(i,...) { x <-  listsimulator(N.trueadmitter=N.trueadmitter,N.falseconfessor=N.falseconfessor,N.withholder=N.withholder,N.innocent=N.innocent,baseline_binomial_prob=baseline_prob)$reject_null}))
    admitter_Nz <- c(admitter_Nz, admitter_Ns[i])
    baseline_probz<- c(baseline_probz, baseline_probs[j])
    powerz <- c(powerz, power)
  }
}

graph.frame.4 <- data.frame(x=baseline_probz, y=admitter_Nz, z = powerz)

powerfig4 <- levelplot(z~x*y, graph.frame.4, cuts=20, col.regions=col.l,
                       colorkey=FALSE,
                       at=depth.breaks,
                       ylab = "Total Number of Admitters (20% False Confessors)",
                       xlab = "Probability p of binomial distribution with 4 trials",
                       main = "Varibility of responses to non-sensitive items",
                       legend = 
                         list(right = 
                                list(fun = draw.colorkey,
                                     args = list(key = list(col = col.l, at = depth.breaks),
                                                 draw = FALSE)))
)

pdf("placebopowerfig4.pdf")
powerfig4
dev.off()


# Figure A1
pdf("powertiles.pdf", height=12, width=12)
print(powerfig1, split=c(1,1,2,2) , more=TRUE ) 
print(powerfig2, split=c(2,1,2,2) , more=TRUE ) 
print(powerfig3, split=c(1,2,2,2) , more=TRUE) 
print(powerfig4, split=c(2,2,2,2) ) 
dev.off()




